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ABSTRACT 

Tidal interaction between a gaseous disk and a massive orbiting perturber is known to result in 
angular momentum exchange between them. Understanding astrophysical manifestations of this cou- 
pling such as gap opening by planets in protoplanetary disks or clearing of gas by binary supermassive 
\ black holes (SMBHs) embedded in accretion disks requires knowledge of the spatial distribution of 

. the torque exerted on the disk by a perturber. Recent hydrodynamical simulations by Dong et al 

(2011) have shown evidence for the tidal torque density produced in a uniform disk to change sign 
at the radial separation of « 3.2 scale heights from the perturber's orbit, in clear conflict with the 
previous studies. To clarify this issue we carry out a linear calculation of the disk-satellite interaction 
putting special emphasis on understanding the behavior of the perturbed fluid variables in physical 
£S| \ space. Using analytical as well as numerical methods we confirm the reality of the negative torque 

density phenomenon and trace its origin to the overlap of Lindblad resonances in the vicinity of the 
perturber's orbit — an effect not accounted for in previous studies. These results suggest that cal- 
culations of the gap and cavity opening in disks by planets and binary SMBHs should rely on more 
realistic torque density prescriptions than the ones used at present. 

Subject headings: accretion, accretion disks — waves — (stars:) planetary systems: protoplanetary 
, disks — solar system: formation 
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Understanding the response of a gaseous disk to a periodic non-axisymmetric gravitational perturbation is an im- 
Ci ', portant astrophysical problem. It arises in a variety of contexts such as the disk-planet interaction in protoplanetary 
disks, orbital evolution of supermassive black hole (SMBH) binaries embedded in circumbinary disks, dynamics of 
accretion disks in cataclismic variables, neutron star and black hole binaries, and so on. 

Gravitational coupling between a gaseous disk and a massive orbiting perturber inevitably leads to the angular 
momentum exchange between the two. This has two important consequences for the evolution of the system. First, 
the orbit of the perturber can change leading to planet migration in protoplanetary disks (Ward 1986) and to the 
black hole inspiral in the case of SMBH binaries embedded in gaseous disks (Ivanov et al. 1999; Gould & Rix 2000). 
Second, the angular momentum lost by the perturber gets absorbed by the disk, which causes redistribution of the disk 
• ■ surface density. In the most dramatic cases, when the perturber is massive enough, gas may be completely depleted 
in some parts of the disk, resulting in gap opening in protoplanetary disks (Ward 1997) and cavity formation around 
^~ 1 ! the SMBH binaries (Armitage & Natarajan 2002). 

It is important to note that the latter process directly affects the former since the strength of tidal coupling between 
the disk and the perturber leading to migration depends on the spatial distribution of the disk surface density (Ward 
1997). Thus, the change of the disk properties caused by the gravitational influence of a perturber has an immediate 
and important effect on the orbital evolution of the perturber itself. 
Disk-plane10 interaction leads to the excitation of nonaxisymmetric density waves, which carry the angular momen- 
\ turn lost by the perturber away from their excitation site. Subsequent dissipation of the waves due to some linear 
or nonlinear process transfers this angular momentum to the disk material, ultimately driving evolution of the disk. 
The torque density (torque per unit radial distance) exerted on the disk material by damping waves dT/dr\d (i.e. the 
amount of angular momentum transferred by dissipating waves to the disk fluid per unit time and per unit radial 
distance) can thus be represented as 
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1. INTRODUCTION. 



dT 
dr 



where dT/dr is the excitation torque density — the amount of angular momentum added to the density wave per unit 
time and per unit radial distance by planetary tide. Operator C describes the wave damping — transfer of the wave 
angular momentum to the disk material by some dissipation mechanism. This operator is intrinsically nonlocal since 
the amount of angular momentum transferred by the wave to the disk at some location depends on the amount of 
angular momentum accumulated by the wave prior to that point. 
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3 From now on we will refer as "disk-planet" or "disk-satellite" to any other type of system where the tidal coupling is present, e.g. a 
SMBH binary surrounded by the disk, an accreting white dwarf-main sequence star binary, etc. 
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Thus, to know how the tidal coupling affects the disk one must understand both the excitation of the density waves 
(i.e. the radial dependence of dT/dr) and the details of their damping (i.e. the explicit form of operator C in equation 
(QJ). This important point was previously highlighted by Lunine & Stevenson (1983), Greenberg (1983), Goldreich 
& Nicholson (1989). The behavior of C has been previously explored for linear viscous dissipation by Takeuchi et al. 
(1996) and for nonlinear damping by Goodman & Rafikov (2001) and Rafikov (2002). The results of these studies 
suggest that for low-mass perturbers wave dissipation is a nonlocal process and most of the wave damping occurs far 
from its driving region. 

In this work we concentrate on the wave excitation and do not consider its propagation and dissipation. Driving 
of the density wave is characterized by the dependence of torque density dT/dr on r, or the radial distance from the 
planet. Spatial distribution of dT/dr has been previously derived from the results of direct numerical simulations of the 
disk-planet interaction (Bate ct al. 2003; D'Angelo & Lubow 2008). At the same time there have been few analytical 
studies of the spatial behavior of the torque density. Linear theory of density wave excitation has predominantly 
concentrated, with few exceptions (Korycansky & Pollack 1993, hereafter KP93; Goodman & Rafikov 2001; Muto & 
Inutsuka 2009), on the behavior of fluid variables in Fourier rather than physical space. 

Nevertheless, some asymptotic results regarding the behavior of dT/dr have been known since the pioneering study 
of the disk-satellite interaction by Goldreich & Tremaine (1980; hereafter GT80). According to this work, far from the 
perturber, at radial separations from its orbit \r — r p \ (r p is the semi-major axis of the circular orbit of the planet) 



Here M p is the planetary mass, £o is the disk surface density assumed to be uniform on scales ~ \r — r p \, fl is the 
angular frequency of the disk at r pi and K n is the modified Bessel function of order n. 

To obtain this expression GT80 have computed the amplitudes of the torque produced by the individual azimuthal 
Fourier harmonics of the planetary potential and then assigned their action in physical space to the positions of the 
corresponding Lindblad resonances. With this prescription the one-sided (i.e. considering only r > r p or r < r p ) torque 
density maintains the same sign independent of the radial separation from the planetary orbit. Torque density behavior 
similar to that given by equation ([2]) was also obtained by Lin & Papaloizou (1979) using impulse approximation. 

The torque density prescription ^ has been extensively used in studies of gap opening by planets (Lin & Papaloizou 
1986; Trilling et al. 1998; Bryden et al. 1999; Armitage et al. 2002; Varniere et al. 2004; Crida et al. 2006) and orbital 
evolution of the SMBH binaries surrounded by gaseous disks (Gould & Rix 2000; Armitage & Natarajan 2002; Lodato 
et al. 2009; Chang et al. 2010; Alexander et al. 2011). Modifications of the prescription ([2]) have sometimes been 
adopted e.g. with the proportionality coefficient different from Cqtso ( e -g- Papaloizou & Lin 1984; Armitage & 
Natarajan 2002). However, the key features of the equation ([2]) namely the \r — r p \~ A dependence and the positive sign 
of the asymptotic torque density for r > r p (and the negative sign for r < r p ) have rarely been questioned (cf. Crida 
et al. 2006). 

Recently Dong et al. (2011) investigated the spatial behavior of the torque density using high-fidelity numerical 
simulations of the disk-planet interaction in the two-dimensional shearing sheet geometry. Quite unexpectedly, they 
found that the one-sided torque density does not maintain a constant sign but rather changes from positive to negative 
(for r > r p ) at \r — r p \ « 3.2/i, in disagreement with the conclusions of GT80. Dong et al. (2011) have shown this 
result to be quite robust and independent of the numerical issues (e.g. resolution of the simulations, simulation box 
size, etc.). 

In this work we demonstrate that some of the linear results of GT80 need to be revised at the qualitative level 
even within the framework of the linear theory. In particular, we show that the torque density does indeed change 
sign and becomes negative for r > r p (positive for r < r p ) beyond several scale heights in radial separation \r — r p \ 
from the planet, in full agreement with the numerical results of Dong et al. (2011). This has important implications 
for understanding the issues of the gap and cavity opening by planets in protoplanetary disks and SMBH binaries in 
circumbinary disks correspondingly. 

Our paper is structured as follows. In <J2] we summarize the linear equations for the perturbed fluid variables, 
which are subsequently solved numerically in <J3] We demonstrate good agreement of the results with the analytical 
calculations based on the global Airy representation developed in £13.21 We compute the spatial behavior of the torque 
density and the angular momentum flux in Sj4] and Sj5] correspondingly . The main result of this paper — confirmation 
of the negative torque density phenomenon — is described in ^4.11 We discuss the origin of the negative torque density, 
performance of the Airy representation, and compare our results with existing studies in |J6J Astrophysical implications 
of our results are outlined in Sj7] 



We study tidal coupling of a planet with a uniform disk in the shearing sheet geometry (Goldreich & Lynden-Bell 
1965), which allows us to neglect geometric curvature effects while preserving the main physical ingredients of the 
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2. BASIC EQUATIONS. 



4 We assume that GT80 calculation refers to dT/dr and not to dT/dr\d since no explicit dissipation mechanism was mentioned in their 
work. 
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system. We also neglect the vertical dimension and assume the disk to be two-dimensional. In this approximation one 
works in a Cartesian coordinate system with x — r — r p and y = r p (tp — tp p ) playing the role of radial and azimuthal 
coordinates correspondingly. 

The dynamics of fluid is then governed by the equations of motion and continuity (Narayan et al. 1987; hereafter 
NGG): 

dv VP 

— + (v ■ V) v + 2VL x v + 4A£lze x = — — - V$, (3) 
BY, 

— + V-(vS)^0, (4) 

where v = (v x ,v y ) is the fluid velocity with components in the x and y directions correspondingly, P is gas pressure, 
$ is the planetary potential, and A = (r/2)(dTl/dr) is the shear rate. When planet is not present and the disk is 
homogeneous, E = So, these equations have a solution in the form of a simple linear shear: v Xt o — 0, v y .o = 2Ax. 

When a planet of mass M p is present (at x = y = 0) a stationary (i.e. d/dt = 0) pattern of the density and 
velocity perturbations gets established in the disk. We can then analyze equations ©-(111) in a standard way, see 
GT80 and NGG for details. First, we assume E = So + Si, v x — i> x ,i, v y = 2Ax + v Vt \ and linearize equations 
for small perturbations Si, v x ,i and V y> i- We then represent all perturbed fluid variables via Fourier integrals as 
{Si, v Xt x, $} = (2tt)^ 1 dk y exp (ik y y) {<5S, u, v, 0}, which takes care of the y-dependence leaving us with a set 
of equations in x coordinate only. 

NGG have shown that the resulting system can be reduced to a single differential equation for the azimuthal velocity 
perturbation v only: 

f AA 2 k 2 , n 2 + klc 2 \ 2Akl 2B 



y -x 2 



a v 
dx 2 

while the perturbed radial velocity u and the surface density perturbation <5S are directly expressed in terms of v via 
the following relations: 

k y c 2 ^- + AABk y xv + 2Bk y ^\ , (6) 



fc2 c 2 + 4B 2 y f da 

5S = ,„ 2 S ° p „ ( 2B^- - 2Ak 2 xv - k 2 A . (7) 
k 2 c 2 + 4B 2 \ dx v v J v ; 

Here c is the sound speed, B = O + A is the Oort's constant, k = (4BQ) 1 / 2 is the epicyclic frequency and <f> is the 
Fourier component of the planetary potential $ given in a point-mass case by 

,,, . GM P 86 , .GM Pl T ^ 

4>{k y ,x) = Ko{\kyx\), — =sgn(x) -k y Ki {\k y x\) . (8) 

7T OX IT 

In a Keplerian disk A = -(3/4)0, B = and k = O. 
Introducing a new dimensionless radial coordinate 



z(x, k y ) = x(4\A\k y /c) 1/2 — ^ y^kyh) ' x(z,k y ) = zh ^-j^-kyh^J , (9) 

where h = c/Q is the scale height, we rewrite equation ([5]) in a standard parabolic cylinder equation form 

d 2 v 



. )j2 +vl--a)=R(k y ,x), (10) 



)= M, _^0 A _ GM^ 



Here 



2c r 2\A\k y cdx 2irc 
_k 2 + k' 2 c 2 _ k 2 1 + (k y h) 2 (n/n) 2 _ 1 + (k y h) 



kyxK (\k y x\) + sgn(a:)— K\ (\k y x\) 



A\A\k y c 4\A\n kyh 3(kyh) 



(11) 



(12) 



where the last equality holds for a Keplerian disk. Clearly, a ^> 1 both when k y h <C 1 and fc y /i 3> 1, i.e. for the modes 
excited at Lindblad resonances lying either very close to the planet, at |x| < h, or far from the planet, at |sc| > h. 
Density wave harmonics with k y h ~ 1, which carry most of the angular momentum flux in a disk of uniform surface 
density and are excited at |a^| ~ h have a ~ 1. 
Equation (JTOj) must be solved with the boundary condition in the form of the outgoing wave for x — ¥ ±oo. 
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3. SOLUTIONS OF FLUID EQUATIONS. 

Solution of equation (jlOp can be directly expressed in terms of the parabolic cylinder functions, but they are not easy 
to analyze. Thus, instead of working directly with the parabolic cylinder functions we solved equation (|10p numerically 
as outlined below in ij3.ll Also, as we describe in £13.21 and Appendix ^ whenever a > 1 one can globally approximate 
the solutions of equation (ITUl) via the combinations of Airy functions of special arguments. We call this approximation 
the Airy representation and it forms the basis for our subsequent analysis of the disk-planet coupling. 

3.1. Numerical procedure. 

Direct numerical integration of the equation ([5]) is based on the procedure developed in KP93, where a shooting 
method is implemented by matching the numerical solution from the origin to the WKB outgoing waves. 

As shown in NGG, the exact homogeneous solution to equation ([5|), i.e. the parabolic cylinder function, can be well 
approximated by the WKB outgoing waves 



v(x^±Oo)~ /_f_ e ±<3*»* */4h (13) 

y 3k y x 

(in a Keplerian disk) when x/h ^> (8/9)(l + (k y h) 2 )(k y h)~ 2 . Additionally, for this approximation to accurately 
represent the inhomogeneous solution of ([5]) one requires \xk y \ > 1, so ^ < 1 and d<p/dx <C 1. Both conditions 
together imply that x/h 3> min{l / \k y h) , 1}, which is the criterion we use to match our numerical solutions to the 
outgoing waves. 

Following KP93, to obtain solutions satisfying boundary conditions we start by shooting two arbitrary linearly 
independent homogeneous solutions and v% and one arbitrary inhomogeneous solution v % from the origin to the 
boundaries at ±xb with Xs/h 3> min{l/(fc a /i), 1}. Then we can write the desired solution as a linear combination of 
these solutions 

v = v* + Oi«i + a 2 v%, (14) 

where ai and a 2 are two complex constants to be determined by the boundary conditions following from the equation 
C3J): 

d 

— {v l (-x B ) + aiv?(-x B ) + a 2 v^(-x B )} 

= yi^kyXs + ^T~ ) {v\-~x B ) + aiVi (-x B ) + a 2 v^ (-x B )} , (15) 
d 

— {v l (x B ) + oiv^xb) + a 2 v 2 l (x B )} 

= (dkyXB - — ] {v l (x B ) + aiV^(x B ) + a 2 v%(x B )} ■ (16) 

V 2 okyX B J 

This method is quite efficient since we only have to integrate three solutions and invert a 2 x 2 matrix to find an 
exact solution for a given k y . Moreover, we have numerically proven that x B /h does not need to be too large compared 
to min{l / \k y h) , 1} to give the right answer with good accuracy. 

In the subsequent numerical calculations we use a 4-th order Runge-Kutta integrator with spatial resolution of h/200 
for 320 uniformly log-spaced values of k y h between 0.01 and 15, and potential softening length of 10~ 4 /i. 

3.2. Analytical Airy representation. 

In Appendix [X] we show that the solution of equation (|10[) with the boundary condition in the form of outgoing 
waves can be expressed in the limit a» 1 via the Airy functions, which have some algebraic function of x as their 
argument, see equations (|A2[) - (|A6|) . There we also derive the expressions (|A13I) and (IA14[) for the real and imaginary 
parts of v, which contain a number of exponentially small terms. Neglecting these terms, which is valid when a» 1, 
i.e. when either k y h 3> 1 or k y h -C 1 one finds that the solution for v reduces to 



TTg(z) 



a 



1/2 



OO Z 

Ai(-t(z)) J g{s)Bi{~t{s))R(s)ds + Bi(-t(z j) J g(s)Ai(-t(s))R(s)d6 



(17) 



ng(z) 



sM—^+^-'W). ( 18 ) 



a 



1+ = / 9{s) 



e 



Ai(-*(*)) - -=-Bi(-f(*)) 



R(s)ds. (19) 



Here Ai and Bi are Airy functions, while functions g and t are defined by equations 
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Fig. 1. — Behavior of the individual Fourier harmonics of the azimuthal velocity perturbation v as a function of both k y (left panels) and 
x (right panels). Respective values of x or k y h are shown in each panel. We plot 5R(i>) and 3(-u) obtained both by the direct numerical 
integration of the equation 101 (solid curves) and using the Airy representation (dashed curves). 

We also demonstrate in Appendix [3] that whenever k y h < 1 and a > 1 the factor I + defined by equation (TIT))) can 
be approximated for a Keplerian disk as 



/+ = R{2V^h(3k y hy 1/2 ) 



l/2\ _ 



GM P 
' 2ttc 



wiTo (M) + (M) 



\Jl + {k y hf 



(20) 



Equations ([T7l) -(|20 |) provide the foundation for our subsequent linear analytical study of the planet-generated density 
waves in §jj3j [5] 



3.3. Results. 

In Figure [T] we show the behavior of the Fourier harmonics of the azimuthal velocity perturbation v as a function 
of x for a fixed k y and also as a function of fcj, for a fixed x. We display the behavior of $t(v) and obtained by 
numerically solving equation (|10[) . and also the analytical Airy representation of the same variable given by equations 

For a fixed k y satisfying the conditions k y h < 1 or k y h > 1 (i.e. for a > 1) the behavior of u is reproduced by 
the Airy representation very well — in most cases the analytical prediction is hardly discernible from the numerical 
calculation, see Figures [TJi,f. But as FigureQJ shows the Airy representation actually works quite well even for k y h ~ 1 
when formally it should not be applicable: analytical prediction falls essentially on top of the numerical results for 
both $t(v) and 9?(u) as long as x > 5/i; significant discrepancy between the two is noticeable only for x < 2h. 

Results for fixed x te nd to confir m this trend. Already at moderately large separations from the planet, e.g. x = 5h, 
analytical predictions (|Al3j) - (|A14j) reproduce the numerical results for all values of k y (including k y h ~ 1) extremely 
well, see Figure [Tfc. Even at x = 2h the agreement is still pretty good, see Figure [TJd, even though at this value of x 
the modes with k y h ~ 1 play the dominant role. It is only very close to the planet, e.g. at x = 0.2h as shown in Figure 
QJi, that the discrepancy between the Airy representation and the direct numerical calculation becomes significant, but 
again only for k y h ~ 1 (analytical formulae agree with the numerical results in the limits k y h <C 1 and k y h 3> 1 even 
at x = 0.2h). 

In Figure[2]we plot the behavior of the surface density perturbation Ei(a;, y) in physical space. We compute Ei fully 
numerically by calculating its Fourier harmonics using equation with numerically determined v and then taking 
the inverse Fourier integr al. W e also calculate Ei semi-analytically by repeating the same procedure with and 
3(u) given by equations (|A13[) - (|A14[) . One can see that the two ways of computing Ei(x, y) agree with each other 
quite well. They do show some discrepancy at small x rs h, which is expected given the dominant role of the modes 
with k y h ~ 1 in shaping the harmonic content of Yi\{x,y) close to the planet. However, further out from the planet 
the agreement between the Airy representation and the numerical calculation improves, see Figure [2]:. Azimuthal 
density wave profiles at different separations from the planet have previously been computed by Goodman & Rafikov 
(2001) with the same method as we employ here, and by Dong et al. (2011) using direct numerical hydrodynamical 
simulations of tidal coupling between the low mass planet and the disk. Their results are in agreement with ours. 

These comparisons clearly illustrate the robustness of the analytical Airy representation for the azimuthal velocity 
perturbation v in certain limits. They also show that this representation can be used for accurate calculation of the 
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y/h- 3/4(x/h) 2 y/h - 3/A(x/h) 2 y/h - 3/A(x/h) 2 



Fig. 2. — Azimuthal cuts through the density wake in coordinate space showing the surface density perturbation Si at different radial 
separations from the planet (x = h,2h,3h as labeled), normalized by x 1 ' 2 . Solid curves show numerical results while the dashed curves 
correspond to Si profiles obtained using Airy representation. An azimuthal shift by (3/4)x 2 is applied in each panel to facilitate comparison. 

behavior of the fluid variables in physical space at large x/h. All this gives us confidence in the subsequent use of this 
analytical approximation for the calculation of the asymptotic behavior of the torque density in 21 and the angular 
momentum flux in SJ5J 

4. TORQUE DENSITY. 

We now proceed to study the spatial behavior of the torque density dT/dx — the amount of torque that is exerted 
by a planet on the disk per unit radial distance x. Analytical calculation of this important physical quantity has been 
first carried out by GT80 but only in the asymptotic regime \x\ > h. Besides, as we show below, this calculation was 
incorrect. Also, a number of numerical studies have derived dT/dx for arbitrary x as a by-product of their simulations 
(Bate et al. 2003; D'Angelo & Lubow 2008). Our current semi-analytical calculation is intended to provide the full 
description of dT/dx in the linear regime for any x. 

The torque density is defined as 



dT 

dx 



ay 



(21) 



Note that compared to dT/dr this definition lacks an extra factor r p and strictly speaking represents the momentum 
density (even though we will still call it torque density in this work). This is because r p is not a well defined variable 
in the shearing sheet setup. 

Fourier decomposing Si and $ in azimuthal direction and manipulating the resultant expression we arrive at the 
following alternative expression: 



dT 
dx 



d l) k dk ^ 



(S)r" 47rWS) 



_87rEn_ 



- Akfx^tv) 
ox y 



(22) 



The last equality was derived using equation (J7J. 

Our numerical solutions for the azimuthal velocity perturbation v obtained in the previous section by numerically 
integrating the linear fluid equation (|10|) allow us to compute the individual Fourier harmonics (dT/dx)k using equation 
(|22p. In Figure[3]we show the behavior of (dT/dx)^ as a function of both k y and x. We also display (dT/dx)k computed 
in the framework of the Airy representation using two levels of accuracy. Dot-dashed curve is derived from equation 
(|2"2")l using v given by equation (|18l) . which neglects the exponenti ally sm all terms, with I + given by equation (fT^jl . 
The (dT/dx)k shown by the dashed line uses v given by equation (|A14I) in which all subdominant terms have been 
retained. 
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Again, we see very good agreement between the theoretical Airy representation and the numerical results in the 
limits of both k y h > 1 and k y h < 1 for all x, as well as in the limit of x > h for any k y h . This is in full agreement 
with the results for v shown in Figure Q] since (dT/dx)k is directly derived from v using ([22]). One can also see that the 
use of equation (|A14[) instead of (1181) to represent v gives a more accurate approximation of the numerical results. For 
that reason we generally recommend using the more accurate equation (|A14|) in semi-analytical calculations involving 
the Airy representation. 

Next we explore the behavior of the full torque density dT/dx, integrated over all k y according to the equation ([22]) . 
In Figure 0] we present the result of such calculation using numerically derived v. We also plot the analytical dT/dx 
computed using Airy representation of v at two levels of accuracy. First, simplified calculation (dot-dashed line) uses 
given by equation (| 18[) with I + computed according to the approximation (|19|) . Not surprisingly this method 
of computing dT/dx does not fare well at x ~ h and it shows significant discrepancy with the numer ical s olution as 
x/h — > 0. Second, more accurate version of the Airy representation (dashed line) relies on equation (|A14j) retaining 
the exponentially small terms (which, however, become important as k y h ~ 1) to describe the behavior of 3(u). This 
calculation yields better agreement with the numerical result as shown in Figure |4] even though there is still significant 
discrepancy (at the level of 10 — 20%) as x ~ h. At the same time, as expected based on the results of fJ3]both versions 
of the Airy representation work quite well at large separations from the planet, at x > 2h. 

We also plot in the Figure [4] dT/dx derived from the direct numerical simulations of the disk-planet interaction 
(Dong et al. 2011). Clearly, the agreement between the linear theory results and the outcome of direct simulations is 
very good. 

4.1. Negative torque density phenomenon. 

Different methods of calculating dT/dx shown in Figure H] universally exhibit the phenomenon of the negative torque 
density first described in Dong et al. (2011). Using direct numerical simulations these authors have shown that in a 
uniform disk dT/dx changes sign at the finite separation from the planet in apparent disagreement with the predictions 
ofGT80. 

To illustrate this behavior in more detail in Figure [5^ we show a zoomed-in view of dT/dx at large separations 
x > 3h. There we plot results of our numerical calculation, analytical calculation using Airy representation with 
3(i>) from equation (|A14I) . and dT/dx derived from the direct numerical simulations of Dong et al. (2001). All three 
methods of computing dT/dx agree on the distance X- ss 3.2h where dT/dx changes sign. Linear theory results (both 
numerical and analytical) show good agreement with each other beyond this point while the dT/dx derived from the 
direct hydrodynamical simulations exhibits some features related to numerical issues (mainly determined by the finite 
extent of the simulation domain in x-direction, see Dong et al. 2011). 

We also show on this plot the asymptotic behavior of dT/dx according to GT80 (thin solid line), see equation @. 
It clearly disagrees with the linear results obtained in this work and with the results of the simulations by Dong et al. 
(2011) both in amplitude and the sign of the effect at x > 1. We provide explanation for this discrepancy in £16.11 

In Appendix [B] we show using Airy representation of v that the asymptotic behavior of the torque density in a 




Numerical 

Eq. (21) with lm(v) from (A14) 
Eq. (21) with Im(i>) from (17)-(18) 
Dong et al. (2011) 



Fig. 4. — Full torque density dT/dx (integrated over all k y ) as a function of x in the vicinity of the planet. Numerical linea r results are 
given by solid line. Results obtained with the Airy rep rese n tati on (equation {22) integrated using S(d) given by equation JA141 1 are shown 
by dashed line, and those using S(u) from equations 118M19| I are shown by dot-dashed line. Dotted curve shows dT/dx extracted from 
one of the numerical simulations of Dong et al. (2011). 

Keplerian disk valid in the limit x/h — > oo is given by 
-J ^sgn(,)C 
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-0.613096. 



(23) 



This expression replaces equation ([2"]l. which previously described the asymptotic behavior of dT/dx according to 
GT80. 

In Figure [SJd we compare the asymptotic scaling given by equation (|23|) with the behavior of dT/dx at large sep- 
arations [x < I2h) computed using our numerical determination of v in equation ()22|1 . We also show on this plot 
dT/dx obtained by numerically integrating the analytical expression (|B4[) over k y . One can see that all three ways of 
representing the behavior of dT/dx at large x agree with each other quite well. 

5. ANGULAR MOMENTUM FLUX. 

To get the full picture of the disk evolution driven by planet-generated density waves it is important to properly 
understand not only the wave damping mechanisms but also the amount of the angular momentum carried by the 
waves — angular momentum flux Fjj, and its spatial pattern. 

Total angular momentum flux carried by the density wave is 



F, 



H 



S / dy v Xt iV Vt 



Fn,kdk y , 



F. 



H,k 



47r£ [3t(v)$t(u) + 3f(u)S(it)] = 47r£ 3f{ (uv* 



(24) 



where Fn.k is the angular momentum flux carried by a particular azimuthal mode of the wave, and * denotes complex 
conjugate. Using equation © and the fact that 4> is purely real we find 



Ffi,k — 



47rSn 



kyC 2 



AB 2 



2Bky<fSf(v) 



(25) 



Because of the angular momentum conservation the radial divergence of the angular momentum flux dFfj ^/dx must be 
equal to (dT/dx)k- To demonstrate this property we differentiate equation (1251) with respect to x and then transform 
the term proportional to v*dv/dx using equation ([S]). In the end we find that dFn,k/dx is indeed given by the same 
expression in terms of ^s(v) as (dT/dx)k, see equation (|221) . 
In the important limit x — > oo we can compute the asymptotic value of the angular momentum flux harmonic Fh,u 
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Fig. 5. — Illustration of the negative torque density phenomenon, (a) Behavior of dT/dx for x > 3h clearly showing change of sign of the 
torque de nsity at x_ « 3.2h. Solid line stands for numerical linear calculation, dashed curve is for Airy representation with Q(v) given by 
equation l|A14| l. and dot-dashed line shows the results of simulations (Dong et al. 2011). Thin solid line shows the asymptotic prediction 
of GT80. (b) Asymptotic b ehav ior of dT/dx at large x/h (note the negative sign). Solid curve is the same as in panel (a), dot-dashed line 
shows asymptotic behavior (123 P . and dashed line results from integrating equation (IB4 1 ) over k y . 



analytically. For this purpose we use equation (|25|) . in which <f> — > as x — > oo (see equation (fTTj) ). so that 



F H ,k{x -t oo) 



4nY, k y c 2 
k 2 y c 2 + 4B 2 



M(v)^(v)-Z(v)-^$t(v) 



(26) 



In the same limit x, z — > oo equation (jTTJ) reduces to 



»(«)~-^/+Bi(-*(2r)), (27) 

where I + is given by equation (|19[) . The exponentially small contributions present in equation (|A13I) have been dropped 
from this expression. 

Using these facts we show in Appendix [C] that the Fourier component of the angular momentum flux is given in the 
limit k y h < 1 and x/h^> 1 by the following formula in a Keplerian disk: 



FH,k{x ->• oo) = -k 



4,_ 2 (GM p ) 2 E 



3 v n 2 



•->A-,. (§)+* (| 



(28) 



This expression coincides with the corresponding result of GT80 obtained in the limit m < r/h (see their equation 
(91)). 

In Figure [5] we present the run of the full angular momentum flux Fh(x) in real space, as a function of the radial 
separation from the planet. One can see that Fh starts with a non-zero value of 0.098 (in units of (G 2 M 2 T,q) / (he 2 )) 
at x — and then gradually rises until x is about 3. Beyond that point Fh{x) starts decreasing with the distance 
although the decline is rather modest: the peak value of Fjj is 0.938 (in the same natural units) while the value of 
Fh at x = 15h is about 0.908. This reduction of Fh(x) with the distance is a direct consequence of the negative 
torque density phenomenon discussed in H4.1[ since dFu/dx = dT/dx as we have shown before. The latter point is 
additionally emphasized in Figure [5] where we also plot the integrated torque T(x) = J Q (dT/dx)dx. One can see that 
Fh(x) and T(x) are identical to each other up to the vertical shift by Fh{0). 

Previously, using calculation in the Fourier space, GT80 have found Fh (oo) = 0.93, which is close to the value we 
find. The small difference between these numbers might possibly be ascribed to the small k y range used by GT80 to 
compute Fh (oo) and should not be taken too seriously. 

6. DISCUSSION. 

Most of the existing linear calculations of the disk-satellite coupling have been limited to the Fourier space (Ward 
1986; Artymowicz 1993; Tanaka et al. 2002). Notable exceptions include Meyer- Vernet & Sicardy (1987) who discussed 
the radial structure of the perturbed velocity eigenfunctions for a fixed azimuthal wavenumber under a variety of 
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Fig. 6. — (a) Scaling of the angular momentum flux Fh with the distance from the planet (solid curve). We also plot integrated torque 
T(x) accumulated by the wave at distance x (dashed curve). Note that Fh(x) has a non-zero value at x = 0, and that dT/dx = dFu/dx. 
(b) C lose -up view at x > 2h to illustrate the slight decay of Fh(x) and T(x) at large £ as a result of the negative torque density phenomenon 

feittSpt'ions about the physics operating in the disk. Later on KP93 have done similar calculation for a disk with 
non-zero pressure (the same setup as used in our current work) and computed the run of Fh,h as a function of x. But 
these studies have only looked at the spatial behavior of individual modes driven by the planetary potential. 

More recently Goodman & Rafikov (2001) have integrated the contributions of separate modes over the azimuthal 
wavenumber to obtain the full spatial distribution of the perturbed surface density. Our current work extends the 
study Goodman & Rafikov (2001) by computing the spatial behavior of other important fluid variables such as the 
torque density and the angular momentum flux. 

The key result of our present study is the confirmation of the negative torque density phenomenon (discovered in 
the numerical study of Dong et al. (2011)) in the framework of the linear theory of disk-satellite interaction, and the 
determination of the asymptotic behavior of dT/dx far from the planet. The discrepancy between our result (|23[) and 
the conventional expression ([2]) not just in amplitude of the effect but also in its sign quite naturally demands an 
explanation, which we provide below. 

6.1. Origin of the negative torque density. 

All calculations of the torque and angular momentum flux excited by the planet in GT80 have been carried out 
in Fourier space for individual harmonics of the planetary potential. These results fully agree with our calculations 
in §S|4] & [5l in particular we successfully reproduce their expression for Fn,k in $5] However, GT80 also tried to 
convert their results on torque exerted by the planet to physical (coordinate) space, and in doing this they had to 
resort to certain assumptions. In particular, they postulated that the action of each potential harmonic, including 
the deposition of a corresponding torque contribution, is localized in the vicinity of a respective Lindblad resonance. 
Under this assumption they computed the torque density at each location in physical space as the simple product 
of the spatial density of Lindblad resonances \dk y /dx\ and the amplitude of the torque component in Fourier space 
Tfc = J^° oo (dT/dx)kdx, corresponding to k y with Lindblad resonance at this particular location (i.e. at k y — 2/(3x)). 
Equation ([2]) represents the result of such a calculation in the limit x/h ^> 1. Not surprisingly, one finds the sign 
of (dT/dx) LR to be the same for any x simply because the Fourier amplitudes of the torque have the same sign 
independent of m — there is no disagreement on this issue between us and GT80. 

What we find in this work is that the Fourier contribution to the torque density of each azimuthal harmonic of 
the potential (dT/dx)k is in fact spread out over a finite portion of the disk, rather than localized in a narrow region 
around each Lindblad resonance. In this situation, the overlap of many Fourier modes at a given point in space is 
inevitable and the calculation of the full torque density in real space represents a rather challenging task, requiring 
direct integration of the equation (1221) . As a result of carrying out this exercise, one finds the distribution of the torque 
density different from ([2]) at large values of x, as clearly demonstrated in Figure [5^. 

5 Analytical calculation of Tfc in the limit k y h < 1 using equation iti-'li and Airy representation l|18j l. 1120(1 for ^(v) results in the expression 
analogous to equation (13) of GT80. 
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This reasoning is best illustrated if we temporarily abandon the shearing sheet approximation and consider the disk- 
satellite interaction in the full cylindrical geometry. Then, as pointed out by GT80, coupling between the planetary 
potential and the disk occurs at discrete Lindblad resonances radially separated from the planetary orbit by XL jm ~ 
2r p /(3m), where m = k y r p in the azimuthal wavenumber of the potential harmonic and an assumption m > 1 has been 
made. The distance between the adjacent m-th and (m+ l)-th resonances is Axl,tu = XL, rn /m = XL,m{h/r p )(k y h)~ 1 . 

We need to compare this distance to the width 5x m of the m-th resonance. According to Airy representation (see £13.21 
and Appendix El the latter is set by the condition t ~ 1 or 2a 2 / 3 (£ — 1) ~ 1 (valid near the Lindblad resonance), see 
equations (| A5|) and (|A15|) . One then finds that the resonance width 5x m is given by Sx m ~ iL. m «~ 2 ' 3 ~ XL^ m {k y h) 2 ^ 
for k y h < 1. Note that 5x m <C XL,m- 

A somewhat more illuminating way of deriving 6x m is based on the dispersion relation for the tightly wound density 
waves (Binney & Tremaine 2008), which is valid for to > 1 or k y h > (h/r): 

^ + k 2 cl = m 2 (n-a p ) 2 ^^-^x 2 , (29) 

where k r is the radial wavenumber of the wave. The exact resonance location x^,m — TpTTi 1 (k/2\A\) corresponds to 
the distance where k r = and the wave couples best to the external potential. As the radial separation from the 
resonance increases k r grows until at \x — XL, m \ ~ Sx m one gets k r Sx rn ~ 1. Beyond this point the wave eigenfunction 
starts oscillating rapidly, and coupling to the external potential becomes less effective. Plugging k r ~ and 
x = XL m + 5x m into the dispersion relation (121))) and expanding it in terms of Sx m <C £L,m one finds Sx m ~ h(kyh)^ 1 / 3 , 
the same expression as that as resulting from the Airy representation. 

Using these estimates of AxL iTn and Sx m one can easily see that Sx m > AxL, m and resonances located at distance 
x from the perturber overlap as long as 



**»(£r ->(£)""■ (30) 

Equivalently, overlap occurs for modes having 



2/5 



3/5 

kyh > ( — ) or m > 



© ' < 31 > 



Interference of different modes takes place even though the width of each resonance Sx m is small compared to the 
resonance separation Xh, m from the perturber's orbit. 

Thus, whenever r p /h 3> 1, which is always the case for protoplanetary and many other kinds of accretion disks, 
resonances overlap in the vicinity of the perturber's orbit (for x satisfying equation (1301) ) and the assumption of 
well-isolated resonances adopted by GT80 for their torque density calculation fails. Given that in our calculation the 
negative values of dT/dx appear only for x > X- s» 3.2h (see i)4.ip one needs the inverse of the disk aspect ratio to 
be r p /h > (x-/h) 5 / 3 7 for the negative dT/dx to appear. This requirement is easily fulfilled in a large variety of 
astrophysical accretion disks, including protoplanetary disks. 

In the framework of the shearing sheet approximation, which corresponds to the limit r p /h — > oo, it is not surprising, 
according to equation (|30[) . that resonances overlap for all x and the isolated resonance approximation adopted in GT80 
never holds. In this limit our asymptotic expression (|23[) works fine for arbitrary x. 

At the same time, far from the perturber's orbit, at separations not obeying the condition (1301) resonances become 
sufficiently spread apart to act in isolation. In this case one would expect the GT80's expression © to be valid. Thus, 
when the full cylindrical geometry of the disk-planet interaction is accounted for dT/dx at fist becomes negative close 
to the planet, for x > ar_, but then should change sign and become positive again (outside of the region defined by 
the inequality ([317|) ) in accordance with the GT80 calculation. The details of this transition can only be elucidated via 
linear calculation of the disk-planet coupling in the cylindrical geometry, which is beyond the scope of this work. 



6.2. Comparison with previous work. 

Quite interestingly, the existence of the negative torque density phenomenon could have been inferred from cal- 
culations done even prior to the work of Dong et al. (2011). In particular, Bate et al. (2003) performed global 
three-dimensional simulations of the disk-planet interaction in cylindrical geometry. As a part of this calculation they 
have computed torque density produced by planets of different masses. Close inspection of their Figure 12 shows that 
at small masses, M p < 0.03 Mj, the torque density changes sign at a finite separation from the planet. This effect 
does not disappear as M p is reduced strongly suggesting that it is a linear phenomenon. 

More recently D'Angelo & Lubow (2008) have computed radial run of dT/dx as one of the by-products of their 
three-dimensional global simulations. In their Figure 7 one can again see the torque density changing sign for low 
enough M p < 10 M®, which is in good agreement with the torque density behavior that could have been inferred from 
the work of Bate et al. (2003). Similar torque density pattern can also be seen in D'Angelo & Lubow (2010). 

Finally, Muto & Inutsuka (2009) have carried out a local linear analysis of interaction between a planet and the 
surrounding viscous disk. Even though their results are affected to some extent by the assumed non-zero value of 
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viscosity (our study assumes viscosity to be zero), their torque density still exhibits negative values beyond x ~ ih 
(see their Figure 5) even at the lowest values of viscosity, in agreement with our work. 

The most likely reason for why the negative torque density phenomenon went unnoticed in previous studies lies in 
the small amplitude of this effect. Indeed, the total negative torque contribution obtained by integrating dT/dx from 
sa 3.2/i to 00 amounts to < 4% of the total integrated torque T(oo) = J °° (dT/dx)dx. Thus, it is not surprising 
that this real effect was previously ignored, most likely on the basis of it being due to some numerical issues (poor 
resolution, small simulation box size, etc.). 

Nevertheless, the apparent presence of the negative torque phenomenon in the aforementioned earlier works provides 
additional support for the robustness of the results presented in this paper. In particular, one may wonder whether 
our conclusions based on a linear calculation performed in two-dimensional shearing sheet geometry would hold when 
some of our assumptions — local approximation, linearity of the perturbation, two-dimensional geometry — are 
relaxed. In this regard we simply point out here that both Bate et al. (2003) and D'Angelo & Lubow (2008) ran three- 
dimensional, global simulations for a variety of planetary masses and one can still discern the presence of the negative 
torque phenomenon in their results. This gives us reassurance in the ubiquity of this rather subtle but important (see 
JF} feature of the disk-planet interaction even in more sophisticated and realistic geometries. 

6.3. Accuracy of the Airy representation. 

Another interesting result of our study concerns the accuracy with which the analytical Airy approximation described 
in £13.21 and Appendix |A1 characterizes the details of the disk-planet interaction. Approximation of the behavior of the 
perturbed fluid variables in terms of Airy functions has been used by many authors before (Artymowicz 1993; Ward 
1986) but only to represent the behavior in the vicinity of Lindblad resonances. As such the Airy representation served 
only as a local analytical tool for computing the azimuthal Fourier components of the density wave properties such as 
the angular momentum flux. 

In our study we use Airy representation in a global sense, to describe the spatial variation of the fluid variables for 
the values of x not limited to the locations of the Lindblad resonances. This procedure works well for describing the 
behavior of the Fourier harmonics of the azimuthal velocity perturbation v and of the torque density dT/dx at any 
x as long as the wavenumber of a particular mode satisfies the condition a > 1, which according to equation (|12p is 
true whenever k y h < 1 or k y h > 1. Moreover, for x > h this approximation is quite reliable even for k y h ~ 1. As a 
result the Airy representation successfully reproduces the asymptotic behavior of dT/dx and Fs,k as a function of x 
in the limit > 1 as we have shown in ij4.ll and <JSJ And as Figures [2] and 0] demonstrate one can also use Airy 
representation to reproduce with decent accuracy the density wave properties in physical space. 

These conclusions are in agreement with the results of Heinemann & Papaloizou (2011) who found that the WKBJ 
approximation valid whenever a > 1 works very well for describing the evolution of the spiral density waves excited 
by the turbulence in the disk. 

7. ASTROPHYSICAL IMPLICATIONS. 

Understanding the negative torque density phenomenon (see ^4. 1|) in uniform disks is the major result of our work. 
It is natural to ask whether this effect can somehow change global consequences of the disk-satellite interaction in real 
disks. 

It is known that planets, which are not massive enough to open a gap around their orbit can migrate radially because 
of the small asymmetry of their gravitational coupling to the inner and outer parts of the disk (GT80); this is the 
so-called Type I migration (Ward 1997). In our shearing-sheet setup such asymmetry is absent by construction and 
planets do not migrate. But the negative torque density is still going to be present in real disks (which are close to 
being uniform locally, near the planet) and one may wonder whether it may affect the speed of Type I migration. The 
answer is no, because the speed of migration is insensitive to the spatial structure of dT/dx and depends only on the 
difference of the full one-sided torques exerted on the inner and outer portions of the disk. As we have shown in f}5]the 
one-sided torque is not changed in our calculation compared to the GT80 prediction, so the speed of Type I migration 
stays the same as long as the disk is roughly uniform. 

Deposition of the angular momentum carried by the density waves drives evolution of the disk away from its initially 
uniform state. This evolution depends on the spatial pattern of dT/dr\d (see equation ([T])) representing the wave 
angular momentum that is being deposited in the disk material. According to equation (fl]) the latter quantity depends 
not only on dT/dx which describes the driving of the waves but also on the process that leads to their damping and 
transfer of the wave angular momentum to the disk material (Takeuchi et al. 1996; Rafikov 2002). However, once the 
latter is understood (i.e. the explicit form of operator C in equation (fTJ) is determined) the disk evolution is going to 
be fully determined by the dependence of dT/dx on x, and thus should be directly affected by the negative torque 
density phenomenon. 

The influence of this effect on dT/dr\d should depend on the details of the damping process. Since the negative 
torque density is a rather small effect (the full amount of negative torque is just several per cent of the full torque, see 
f}5]) it should not affect dT / dr\d if the damping is weak and the characteristic distance over which wave is dissipated is 
Id 3> h. What happens in the case of strong damping, when Id becomes comparable to the distance between individual 
Lindblad resonances, is not so obvious. In particular, it is not clear whether dT/dr\d should become negative at 
|i| > X- as does dT/dr, or whether dT/dr\d will tend to converge to the GT80 prediction ([2]). This issue can be 
addressed in the future by calculations making explicit assumptions about the damping mechanism. And one should 
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keep in mind that strong damping can affect not only propagation of the waves but also their excitation in a non-trivial 
way (Muto & Inutsuka 2009). 

It is worth mentioning that, with the rare exceptions of the works by Takeuchi et al. (1996) and Rahkov (2002), 
in the majority of studies of the perturber-driven disk evolution the wave damping process is essentially ignored and 
it is assumed that the amount of angular momentum transferred to the disk material is given simply by dT/dr (e.g. 
Lin & Papaloizou 1986; Chang et al. 2010), i.e. the assumption dT/dr\d — dT/dr is usually made. However, the 
circumstances under which such assumption is justified have never been explicitly clarified^. In general the correct 
approach to understanding the disk evolution due to tidal interaction must incorporate both the correct form of dT/dr 
(or dT/dx) and the description of the wave damping. If the latter is ignored then according to Goldreich & Nicholson 
(1989) the disk should not evolve at all. 

Knowledge of the spatial behavior of the torque density dT/dx (or dT/dr) produced by a perturber is one of the 
key ingredients for understanding the process of gap opening. In particular, in steady state the shape of the gap is 
uniquely determined by the balance of the viscous torque density in the disk and the radial divergence of the angular 
momentum flux carried by the density waves dT/dr\d- 

A proper calculation of the gap shape is important not only in itself but also for computing the full torque exerted 
by the perturber on the disk. If this torque is non-zero (e.g. due to the disk being absent on one side) then the 
perturber will migrate (Type II migration in presence of a gap, see Ward (1997)) as a result of the angular momentum 
conservation. The speed of migration depends on the shape of the gap since the surface density profile is one of the 
ingredients determining how much net torque the perturber deposits in the disk (Ward 1997). 

Calculation of the excitation torque density dT/dx in a non-uniform disk must be modified compared to that in the 
constant S disk. This is usually accomplished by computing dT/dx\ nu in a non-uniform disk according to the following 
prescription: 



dT(x) 



dx 



Yj(x)- ^ ^ 



Eo dx 



(32) 



where dT/dx\ u is the torque density in a uniform disk, a quantity which is the subject of the calculations done in 
GT80 and in our work (where we call it just dT/dx). One normally uses the GT80's prescription ([2]) for dT/dx\ u (Lin 
& Papaloizou 1986; Trilling et al. 1998; Armitage & Natarajan 2002; Armitage et al. 2002; Lodato et al. 2009) or a 
modified version thereof. Given that in this case the torque density has a definite sign, the tidal interaction between 
the perturber and the gap edges is repulsive and the planet drives a positive angular momentum flux through the disk, 
which is consistent with the gap opening picture. 

Our revision of the dT/dx calculation and the discovery of the negative torque density phenomenon might be thought 
of as giving rise to an interesting paradox. With our new expression (I23|) for the torque density one may think that 
equation (f3"2"|) predicts dT/dx\ nu < for x > x_ « 3.2ft,. Then the application of the prescription (|3"2")l for computing the 
full torque in a broad gap with width exceeding 2x_ would result in the negative angular momentum flux accumulated 
by the wave and the attractive interaction between the planet and the fluid at the gap edges. Apparently, under such 
circumstances the gap would not be able to appear in the first place. 

The resolution of this apparent paradox lies in the use of the prescription (|32p to describe the torque density in an 
inhomogeneous disk. Indeed, the calculation of dT/dx\ u is based on the solutions of the perturbed fluid equations (0- 
([7]) derived under the assumption of a uniform density. However, in a non-uniform disk these equations get modified 
and new terms caused by the gradients of E appear in them. These terms are especially important at the steep edges 
of gaps or cavities. Their presence significantly modifies the structure of the velocity eigenfunctions and makes the 
prescription (|32[) irrelevant. In general one should not expect equation (|32p to provide correct results even for the 
excitation torque density dT/dx (not mentioning the torque density dT/dx\d deposited in the disk by the density 
wave) . 

A fully self-consistent calculation of dT/dx\ nu should be based on equations explicitly incorporating disk non- 
uniformity and thus has little to do with dT/dx\ n . An example of such a self-consistent linear calculation of dT/dx\ nu 
is provided in Petrovich & Rafikov (2012; in preparation) where it is shown that a proper account of the disk non- 
uniformity in fluid equations results in the positive angular momentum flux produced by the perturber and the repulsion 
of fluid away from its orbit for gaps of arbitrary width. 

Knowledge of the behavior of the torque density dT/dr is also important for determining the surface density profile 
in cavities or gaps formed in accretion disks by embedded binary SMBHs. In the case of extreme mass ratio inspirals 
the gap formed by the lower mass BH can be quite narrow (Chang et al. 2010) so that the results obtained in the local 
shearing sheet approximation remain valid. In particular, the negative torque density phenomenon may be relevant 
for this type of objects if the mass of the secondary is insufficient to significantly perturb the disk surface density and 
open a gap. 

However, in the case of SMBHs of comparable mass the binary usually clears out a large cavity with inner radius 
comparable to the binary separation in the disk around itself (MacFadyen & Milosavljevic 2008; Shi et al. 2011). In 
that case understanding the torque density distribution requires one to carry out a global calculation of the density 
wave excitation by the binary in a highly non-uniform disk in full cylindrical geometry. But even in this much more 
complicated setting one may still use the intuition gained in our present study. In particular, equation (|3"0)l suggests 

6 It is clear that this assumption is not appropriate in the case of weak wave dissipation (Goodman & Rafikov 2001; Rafikov 2002). 
Whether it may be reasonable in the case of strong wave damping remains to be explored. 
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that the negative torque phenomenon is not going to be important since the range of x where the resonances overlap 
falls inside the cavity. At the inner edge of the cavity, at separations of order the binary semi-major axis the resonances 
are well isolated and procedure adopted in GT80 to compute dT/dr may actually work. 

On the other hand, if the cavity is very wideQ then the torque is generated by only a handful of individual modes 
since only a small number of low-m Lindblad resonances are going to be located in the region with appreciable gas 
surface density. As a result, the pattern of dT/dx may be more reminiscent of the (dT/dx)k corresponding to a single 
azimuthal Fourier harmonic of the potential, which usually exhibits complicated oscillatory pattern, see upper right 
panel in Figure |3] and KP93. This expectation agrees well with the results of numerical simulations of MacFadyen & 
Milosavljevic (2008), Cuadra et al. (2009), Farris et al. (2011), and Shi et al. (2011) who have found oscillating dT/dx 
in disks with large cavities cleared by the SMBH binaries. 

8. SUMMARY. 

We have carried out linear calculation of the gravitational coupling between the uniform gaseous disk and an 
embedded massive perturber. Our work goes beyond similar earlier studies by putting emphasis on understanding 
the behavior of the perturbed fluid variables in real space as opposed to Fourier space. This allows us to address the 
nature of the recently discovered negative torque phenomenon and explain its origin fully in the framework of linear 
theory. 

We come up with the global analytical representation of the spatial variation of the perturbed fluid variables in terms 
of the Airy functions, which is shown to work quite well for describing the properties of the planet-generated density 
waves both in Fourier and in real space. Using this approximation we derive the asymptotic behavior of the torque 
density far from the perturber and show that it matches the results of numerical calculations. Our linear calculations 
confirm that the torque density changes sign at a finite separation from the perturber in full agreement with the results 
of the direct hydrodynamical simulations. This provides an important correction to the calculations of GT80. 

These results have broad implications for understanding tidal interaction of planets with circumstellar disks and of 
binary SMBHs with the circumbinary accretion disks. In particular, our calculations provide a pathway for deriving 
the spatial torque distributions in such disks, which is important for understanding the issues of gap opening by planets 
and gas clearing by binary SMBHs. 
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APPENDIX 
AIRY REPRESENTATION. 



The following results extensively use the relations presented in Chapter 19 of Abramowitz & Stegun (1972; hereafter 
AS). Homogeneous parabolic cylinder equation 



d 2 v 
dX 2 



x- 



0. 



(Al) 



i.e. the equation (fTU|) with R set to zero has two independent homogeneous solutions, which behave as purely outgoing 
waves at infinity (Narayan et al. 1987): 



v + = E(a, X) = k- 1/2 W(a, X) + ik 1/2 W(a, —X), 
u_ = E*(a, -A) = k- 1/2 W(a, -A) - ik 1/2 W(a, X), 
k= y/l + e 2 ™ - e™, 



(A2) 
(A3) 

where W(a, X) is the standard solution of equation (|Al|) . see chapter 19.17 of AS. One can easily show that the 
Wronskian A = W {E(a, A), E*(a, -X)} = 2e™ is independent of A. 

In chapter 19.20 of AS it is shown that when a is large and positive (a ^> 1) the parabolic cylinder function W(a, X) 
can be expressed via Airy functions for any X > in the following way: 



W(a, X) ~ V^y^l 5 (X)Bi(-i), W(a, -X) ~ 2^t£^ <?(A)Ai(-i), 



where 



r(0- 



X 



t(0 = (4a) 2 / 3 r(£), g(X) 



2/3 



arccos 



2/3 



1/4 



(A4) 



(A5) 



Q\ 2 / 3 / , v 2/3 

g) (ev/^T-arccoshe) , £ > 1. 



(A6) 



We will also need to know the asymptotic behavior of Airy functions as A — > oo (here C = (2/3) A 3 / 2 ): 



Ai(A) . 
Ai(-A) 



2v^AV4' 
cos(C — 7r/4) 

V^aV4 



Bi(A) 



V^aV4' 



Bi(-A) 



sin(C - 7r/4) 



(A7) 



Full solution of inhomogeneous equation f| 10[) behaving as purely outgoing wave for x — > ±oo is expressed via the 
homogeneous solutions v± as 



z oc 

v(k y ,x) = —-^- J v-(s)R(s)ds — — J v+(s)R(s)ds. 



(A8) 



Here it is understood that z = z(x) and R{s) = R(x(s)), where the corresponding dependencies are given by equations 
® and (HU). 

Now we consider the behavior of the solutions for positive x (and z). For a > 1 and z > one can write using 
equations (fA"2j) -(fA"4 f 



-Tva/2 



• V2e™/ 2 W(a, z) + i±—-W(a, -z) ~ ^(*) pi(-t(z)) + iAi(-t(«))] , 

p-TTo/2 

, ^e™/ 2 ^, -z) - i ^VF(a, z) 

v2 



2e 7m Ai(-i(z)) - »-— _Bi(-t(z)) 



(A9) 



(A10) 



where £ = z/(2y/a) and t(z) implies f calculated with this value of £. 
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Integrals in equation (|A8I) involve integration over both positive and negative z (or x) , so we need the Airy function 
approximation for v± for negative arguments as well: 



v+~V2e™/ 2 W(a,-\z\)+i- 



-Tta/2 



-W(a,\z\) 



—M- z ) 



2e 7ra Ai(-t(-z)) + i——m{-t(-z)) 



-ira/2 

V2 



-W(a,-\z\) 



~ ^$(-*0 - iM(-t(-z))] , 

where now z < 0, £ = —z/( 2y/a), and i(— 2) means i computed with this value of £. 
Plugging the expressions (|All|) - (IA10|) into equation (| A8|) we find that for a ^ 1 and z > 



M(-t(z)) 



g(s)Bi(-t(s))R(s)ds ■ 



g(s)Bi(-t(s))R(s)ds 



9f(w) 



4 

7T0(z) 



g(s)Ai(-t(s))i?(s)ds 



+ Bi(-t(*)) 



g(s)Ai(-t(s))i?(s)ds 



ff(s)Ai(-i(s))i2(s)ds 



g(s)Bi(-t(s)).R(s)ds 



o 1 ^ 
x / 



Ai(-f(*)) 



M(-t(s)) 



-Bi(-i(*)) 



-Bi(-i(s)) 



R{s)ds. 



(All) 



(A12) 



(A13) 



(A14) 



In deriving these expressions we have used the fact that 

/ r?(s) {Ai, Bi} (-t(s))R(s)ds = - / <?(-s) {Ai, Bi} (-t(-s))R(s)ds 

JO J —oo 

because R(— z) = —R{z). 

These formulae contain a number of terms which are exponentially small in the limit a» 1. Indeed, using the 
relations (IA7|) and definitions (| A5|) - (| A6|) one can easily show that even for z = £ = when ( = ira/2 is largest and the 
term oc Bi is exponentially large (and the term oc Ai is exponentially small), terms like e _7ra Bi are still subdominant 
as O (e~ 7m / 2 ). All other terms containing exponential factors are even more subdominant. Thus, we can safely neglect 
such terms for a ^ 1 and z > to get an approximate form (fTT)) - (flT))) of the solution for v. 

We now evaluate the integral /+ appearing in equation (|1T[) for a Keplerian disk in the limit k y h <C 1 when 
a w (Skyh)^ 1 ^s> 1, with R(s) — R {h^kyh)^ 1 / 2 s) . The second term in this integral is exponentially small when 
a ^> 1 and can be dropped for simplicity. The asymptotic behavior (|A7j) of Ai suggests that the rest of the integral is 
dominated by s suc h tha t t(s) ~ 1, since for other s the integrand is either exponentially small or rapidly oscillates. 
From the definition (|A5|) of t one can see that this requires t«1, resulting in integral being dominated by the small 
vicinity (|£ — 1| <C 1) of s = 2^/a s» 2(3k y h)~ 1 / 2 > 1. One can easily show that in this region 



g-1 
10 



(A15) 



Using the fact that Ai(s)ds = 1 and that the driving term R(s) varies only weakly in the region where | — 1 1 -C 1 
we find for x < (retaining in (|A15[) only the terms of the lowest order in £ — 1) that I + is given by equation (|20l) . 



TORQUE DENSITY. 

We now proceed to calculate the torque density dT/dx. Using equations (fTE)) , (flT)]) . & (|22l) we write 



dx~) k = k 2 y c 2 + 452^^72 



2? A [s(*)Ai(-t(z))] - Afc 2 ^(z)Ai(-i(z)) 



(BI) 
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This expression is accurate up to the exponentially small terms. 

To perform differentiation w.r.t. x in this equation we use the following relations, which can be easily derived from 
equations (|X5 )) -(|X6 ]) : 



dr _ 1 /£ 2 - 1 dt(z) 
dg(z) 



-(Zkyh) 1 ' 2 



t dx 

V 



'a 



(3k y h)V 2 g(z) 



dx 1 h 

Then we find for a Keplerian disk 

'dT\ _ 2ir 2 Y, n 2 k y I + <t> 
~dx J k ~ k 2 c 2 + n 2 /4: aM 2 c 



4(e 2 - i) 

3(k v h)(k y x)g(z)Ai(-t(z)) 



(B2) 
(B3) 



y/1 + {k y hf 
9{z) 



(g(z)) 2 {j-a(g(z))- 
4a £ 2 - 1 



-Ai(-i(z)) + Ai'(-f(«)) 



(B4) 



where Ai'(X) = dAi(X)/dX. This equation is valid for any x (or z) as long as a > 1. 

To get the full torque density one needs to substitute {dT /dx)k into equation (|21[) and perform integration over k y . 
In general one cannot expect the expression (|B4|) to represent (dT/dx)k over the full range of k y since it is valid only 
for a > 1 and formally fails for k y h ^ 1. 

However, this is not a problem for obtaining the asymptotic expression for the torque density dT/dx far from the 
planet, at x/h 3> 1, which we do next. Indeed, {dT /dx)k is always proportional to the planetary potential <f>, which 
according to equation ([8]) behaves as 4> oc exp [— (k y h)(x/h)] <C 1 for x/h ^ 1 and fcj,/i ~ 1 (this follows from the 
asymptotic behavior of the modified Bessel functions at large values of the argument). As a result, the contribution 
of harmonics with k y h ~ 1, for which a ~ 1 and analytical expression (IB4[) fails, is exponentially suppressed and we 
need not worry about them. For the sam e reason one can safely neglect the contribution to dT/dx produced by modes 
with k y h>l even though equation (|B4I) works well for them. 

We thus expect that for x/h > 1 torque density dT/dx is contributed mainly by modes with k y h < 1 for which is 
not exponentially suppressed, a sa (ikyh)^ 1 > 1 and the expression (|B4[) can be used (this also allows us to neglect 
a number of terms 0(k 2 h 2 ) in (|B4[) ) . This expression contains terms proportional to either Ai or Ai' and we look at 
their corresponding integrals (dT/dx)1 s and (dT/dx)^ s separately, starting with the former. 

Airy function Ai(X) rapidly oscillates for X — > — oo and decays exponentially for X — > oo, see equation (| A7|) . In both 
limits the corresponding contributions to the integral are small and the main contribution is provided by the interval 
of k y in which X = —t(z) ~ 1. From equation (IA5|) and the fact that a > 1 it follows that t ~ 1 requires r(£) <C 1, and 
|£ - 1| < 1. Then according to equation (|ATo|) -t(z) 2a 2/3 (l - = a 2/3 (2 - ik y x) and it is clear that {dT/dx)l 8 
is dominated by a narrow vicinity of k y = (2/3)x _1 (with absolute width of 5k y ~ |a;| _1 a~ 2 / 3 ~ |a;|~ 1 (/i/|a;|) 2 ' /3 and 
relative width Sk y /k y ~ (/i/|x|) 2 / 3 <C 1). Inside this narrow range of fcj, all other factors in the integrand function vary 
only weakly and can be evaluated for k y (2/3)x~ 
one finds that 



dT\ as 
dx J 1 



296 [2 

An - 

405 V 3 



and g given by (|A15|) . and taken out of the integral. As a result 
(GM P ) 2 £ 1 



2Ka - 



K, - 



n 2 



(B5) 



where equations ([5]) and (12^1) were used to evaluate <j> and I+- 

Evaluation of the torque contribution (dT/dx)™ proportional to Ai' is a more intricate exercise. It is accomplished 
via integration by parts, which again results in the integral of a rather complicated function multiplied by Ai(—t(z)). 



After laborious but straightforward calculation this function can be evaluated at k y 
becomes 



(2/3)x 1 and the whole integral 



dT 

dx 



16 
27 



2 ± K (~\ Ik (- 

10 °V3/ 3 1 3 



2K n - 



k ; 



(GM p ) 2 So 1 

n 2 x i 



(B6) 



Summing up (IB5|) and (|B6[) we arrive at the final asymptotic expression (|23p for the torque density valid in the limit 
x/h > 1. 

ANGULAR MOMENTUM FLUX. 
Plugging in (fl8]> . ([T9"} and ((23) into ((26]) we obtain 



4irY, k y c 2 
k 2 y c 2 + AB 2 

n 2 P + g(-z) 



Ai(-t(*))^ [ff(*)Bi(-t(z))] - Bi(-i(z))^ [ 5 (z)Ai(-t(z))] 



(CI) 
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As x — > oo so does 2 — ► oo for any Then, according to equations (|A5|) - (| A6|) one has £ — > oo, t(£) — > (3£ 2 /8) 2 / 3 , 
— » (3 a£ 2 / 2) 2 / 3 , and g — > (Sa/2^) 1 / 6 . Also, in equation (|C1[) we can replace Airy f unct ions with their asymptotic 
behavior (|A7[) in terms of trigonometric functions. Plugging all that into the equation (|C1[) . taking the limit fc y /i < 1, 
and using the expression (|20|) for / + one finds for a Keplerian disk 



F H ,k(x -> oo) 



12fe 



■ (GM P ) 2 S 

ft 2 



1 



«>#b(H) + ,-MM) 



(C2) 



with w defined in equation 
equation ([28| . 



20]) . In the long wavelength limit k v h < 1 one has w — ?> 2/3 and equation (|C2j) reduces to 



